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ABSTRACT 

Accurate stress and strain calculations at a notch usually require a non-linear finite 
element analysis when local yielding has occurred. The strain energy density hypothesis 
is a method to predict these stress and strain values. This method proposes that the 
plastic strain energy density is equivalent to the strain energy density found assuming the 
material to be entirely elastic. This hypothesis was evaluated using the finite element 
method, which was tested by comparing to exact solutions of elastic and elasto-plastic 
problems, to calculate the stress and strain field for two notched plates of varying widths 
under elasto-plastic loading. For both geometries, a plane stress and plane strain analysis 
was performed. 

The elasto-plastic strain energy density from the finite element method was found 
to be greater than that predicted by this proposal, which in turn resulted in under- 
predicting the local stresses and strains. This difference was greater for the plane stress 
condition than for the plane strain condition. Comparisons were also made with notch 
stresses based the Neuber method. The two methods appear to give an upper and lower 
bound to the actual stresses and strains. By combining the results of the strain energy 
density method and the Neuber method, reasonably accurate estimates of stress and strain 
values can be obtained. 
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I. INTRODUCTION 



A. FATIGUE LIFE 

Fatigue phenomena were first considered in 1850 when it was discovered that 
railroad car axles failed after a certain amount of time under cyclic loads. The concept of 
endurance limit was introduced, which states that below a minimum load cyclic load, or 
endurance limit S e , some materials will never fail [Ref. 1]. To prevent fatigue failure, one 
only has to design the structure strong enough so that the stresses remain below this 
endurance limit. However, certain structural applications require low weight, and this 
excessive design would be unfeasible. Such is the case with aircraft. All aircraft are put 
under cyclic loads every flight - from repetitive takeoffs and landings, basic 
maneuvering, and flight gust loads. While aircraft are designed to withstand mechanical 
failure modes of large scale yielding and sudden monotonic fracture, it is unrealistic to 
design an aircraft that will never experience fatigue. 

It is not surprising then that the service lives of U. S. Naval Aircraft are based on 
the fatigue life of the aircraft structure. Due to both a safety factor and a cost factor, the 
structural life of every U. S. Naval aircraft is thoroughly tracked. Engine or aircraft 
system components can be replaced, but when the main structural components of an 
aircraft fail, there are few options. The aircraft is either taken out of service, or major and 
costly rework such as replacing the wings on the A-6E Intruder aircraft is required to 
extend its life. Known weak spots are constantly inspected and reworked, if required, to 
extend the life of the aircraft. These critical areas of failure, or crack initiation, usually 
occur at stress concentrations, such as those occurring at notches or rivet holes. Current 
fatigue calculations can estimate the cyclic life of a component given the range of stresses 
and strains it will undergo. However, for these calculations to be accurate, precise stress 
and strain figures are required. This is especially true when the cyclic fatigue calculations 
occur over tens of thousands of cycles, and any error gets multiplied many-fold. Thus, 
accurate calculations of the stress at a notch root is the first step in accurate fatigue life 
estimation. 
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B. STRESS CONCENTRATION FACTOR 



Stress concentrations are stresses that are locally higher than the far-field stress for 
any loaded material. Stresses increase and decrease due to the geometry of a part, and 
when abrupt geometric changes occur, such as a notch or a hole, stresses can be two to 
three times greater than the far-field stresses. Given a defined geometry and applied 
loads, a ratio of local (notch root) stress a to the far-field stress S can be determined and 
remains constant as long as a and S remain in the linear range of the stress-strain curve. 
For a notched specimen, the maximum stress will occur at the notch root, and the ratio of 
the maximum stress to the nominal stress is the stress concentration factor: 




Throughout this thesis, K, is the theoretical stress concentration factor based on a 
elastic model. Even though it is labeled theoretical, its value may be obtained from 
analytic solutions or a finite element analysis. Regardless, it is a constant and holds for 
any stresses in the elastic range. 

The difficulty in calculating the stresses for a geometry such as a notch or a hole 
arise when local yielding occurs. Although an analytic solution or a finite element 
analysis can easily calculate the stress concentration factor that remains constant in the 
elastic range, as soon as plastic yielding occurs, this stress concentration factor decreases 
in magnitude as the yield zone around the notch increases. This may often be the case, 
since local yielding is allowed in the design of the aircraft, and for stress concentrations 
of two to three, it is difficult to design a structurally efficient wing that will prevent local 
yielding. Once the local stress departs the linear range of the stress-strain curve, it thus 
becomes more difficult to calculate the local stresses and strains at the notch root. A 
finite element analysis is possible, but this becomes expensive and time consuming when 
used in high-cycle fatigue calculations that involve tens of thousands of runs. 

In 1961, H. Neuber [Ref. 2] derived a relationship for determining stresses and 
strains at a notch that has been loaded into the plastic range. Although Neuber’s 
derivation involved a notch loaded in antiplane shear, it has been widely applied to 
general notch problems. By the 1980’s, the Neuber method had been adopted by virtually 
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all fatigue analysts, including the Naval Air Systems Command. However, over the years 
since Neuber’s original paper and ever since it has come into prominent use, many 
investigators have proposed alternative means. 

One proposal by Glinka et al. [Ref. 3, 4, 5, 6] is based upon the concept that the 
strain energy density of the material in the yielded zone is virtually the same as the strain 
energy density considering the material to be elastic. This is represented below in Figure 
1.1, where W e is the strain energy density assuming an elastic material and equals the area 
under the linear curve, and W p is the strain energy density for an elasto-plastic material 
and equals the area under the nonlinear curve. This conjecture results in being able to 
calculate the stress concentration factor in the elasto-plastic zone from the strain energy 
density of the elastic model. This proposal is based on reasoning that the for local plastic 
yielding, there is a relatively large volume of material in the elastic region surrounding 
the plastic zone. Glinka continued to work on his proposal through the 1980’s, and 
published several papers applying his model to plane stress and plane strain problems. In 



a 




Figure 1.1. Representation of Strain Energy Density Equivelance Concept 

1992 W. N. Sharpe, Jr., C. H. Yang, and R. L. Trengoning [Ref. 7] evaluated the Glinka 
relations with that of Neuber’s for various plane strain and plain stress configurations 
using experimental data at the notch root. Their conclusions were mixed, and stated that 
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some cases were better predicted by the Glinka model, while others were accurately 
predicted by the Neuber model. 

Drawing upon advanced finite element techniques, the proposal that the strain 
energy density in the plastic zone is equal to that calculated on the basis of an elastic 
solution has been further tested, not only at the notch root, but throughout the plastic zone 
of the model. Elastic and Elasto-plastic problems with closed form solutions, along with 
previously published experimental data have verified the finite element modeling, which 
was then used to calculate strain energy density. Assessments are then made of its impact 
upon fatigue life calculations of aircraft as compared to the Neuber approach. 
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II. THEORETICAL BACKGROUND 



This thesis will concentrate on the Glinka model for calculating the stress 
concentration in the elasto-plastic range. The Neuber model is reviewed below, and will 
be used for comparison of results at the notch root with that of the Glinka model. 



A. 



NEUBER’S MODEL 



Neuber proposed that at the notch root under plastic yielding, the elastic stress 
concentration factor was the geometric mean of the stress concentration factor and strain 
concentration factor, as shown in Equation 2.1 [Ref. 2], 

K,=JkX ( 2 . 1 ) 

where K n = — and K= — 

° S E e 

When the far-field stress S is in the linear rage, this can be rewritten as: 

K] = K 9 K t 

*, 2 = °- 
S e 

Kj^ 



S S 



(2.2) 



or: (K,S) 2 =Eoe 

There are two unknowns here, a and e. To solve for them, the stress-strain constitutive 
relationship is required. One of the most common stress-strain relationships that can be 
applied is the Ramberg-Osgood curve. For a uniaxial stress state, this is: 



a 

£ = — + 

E 



\K 



(2.3) 



where E, K, and n are obtained from a curve fit to the uniaxial stress-strain curve. 
Substituting Equation 2.3 into Equation 2.2, the Neuber relation gives Equation 2.4 
shown below. 



{K,S) 2 _ a 2 



, 



+ a 



K 



(2.4) 
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For the above derivation, the stress-strain relationship using the Ramberg-Osgood 
relationship is based on true stress, a , and true strain, e . However, the computational 
finite element analysis is based on engineering stresses and strains, a and e. The 
relationship between true and engineering stresses and strains, valid up to necking, is 
given below in Equation 2.5. It can be seen that for small strains such as e < 0.01, the 
difference between true and engineering strains is less than 0.5%, and the difference 
between true and engineering stresses is less than or equal to 1.0%. Therefore, for small 
strains, as is the case in this thesis, the true and engineering stresses and strains can be 
interchanged without any significant error. 

£ = ln(l + £) a = a(l + £) (2.5) 

While Neuber’s rule has been well established as an engineering tool to calculate 
notch stresses and strains, it has been shown to overestimate these values [Ref. 3], The 
accuracy of strain estimation is critical for fatigue calculations, therefore, a more accurate 
method would prove very beneficial. 

B. GLINKA MODEL: STRAIN ENERGY DENSITY APPROACH 

Glinka proposed that the energy density at the notch root, calculated on the basis 
of elasto-plastic constitutive laws, is equal to that based on linear elastic constitutive laws 
for equivalent far-field loading. In Chapter III, a detailed derivation of the strain energy 
density will be shown, for now it is given below for the uniaxial case: 

W 0 =^<5dz ( 2 . 6 ) 

Using the above definition of strain energy density, the strain energy at the notch and far- 
field regions are calculated (using a linear elastic stress-strain relationship a = Ee ). For 
the notch root, this becomes: 



W„' £e<fe 


(2.7) 


II 


(2.8) 
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(2.9) 



W a = — 
° 2 E 



For the far- field region, a = S and e = e: 



= E — 

s 2 



( 2 . 10 ) 



Wc = 



2E 



( 2.11 ) 



Rewriting the stresses in terms of strain energy density and substituting these values into 
the theoretical stress concentration factor in Equation 1.1, one gets: 



K = 



K 

w c 



( 2 . 12 ) 



However, Glinka’s hypothesis is that the strain energy density at the root will result in the 
same value, regardless if calculated for a linear elastic or an elasto-plastic material. 
Therefore, this ratio remains constant, even when local yielding occurs at the notch root. 
The argument is that if the area of local yielding is small, and is surrounded by a large 
volume of elastic material, then the energy distribution does not change significantly, 
even when local yielding occurs. 

For a nonlinear stress-strain relationship, W a is found by manipulating the 
integrand of Equation 2.6: 

ad£ = edc + Gd£-edc = d(ce)-edc (2.13) 

Substituting Equation 2.13 into Equation 2.6 results in: 

W a = J[rf(ae)-£<ia] = ge- J q £ da (2.14) 

Substituting the Ramberg-Osgood stress-strain relationship (Equation 2.3) into 
Equation 2. 14 results in the strain energy density in terms of the local uniaxial stress: 



W„ = — + ' 



G0W'-j 



r G_ 



dG 



W= — + 






— + [ — ] (o) ; 
2 E l K 1 V ' 



\n +l . 



(2.15) 



(2.16) 
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(2.17) 



Substituting 2.17 back into the ratio for K { (Equation 2.12) one obtains an equation in 
terms of the far-field stress S, material properties E, n, and AT, and the theoretical stress 
concentration factor K t that can be solve numerically for the local stress a: 



This expression applies only for the uniaxial case. In Chapter III, the strain energy 
density is expressed in terms of the general stress-strain equations. 

By inspecting Equations 2.18 from Glinka and Equation 2.4 from Neuber, one can 
observe that the only difference is the factor of 2/(l+n) in the strain energy density 
model. Since n is less than 1, this term is greater than 1, and for the left side of these two 
equations to be equal, the local stress in Equation 2.18 must be less than the local stress in 
Equation 2.4. Likewise, if the local stress in the strain energy model are lower than those 
in the Neuber model, so will be the local strains. In fact, Glinka states that the Neuber 
model has been shown to overestimate the local stresses and strains [Ref. 3]. 



(K,sy _ a 2 | 2a ( oY 



(2.18) 



E E l+nyK) 
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III. STRAIN ENERGY 



A. STRAIN ENERGY PRINCIPLE 

For a general deformable body, external applied loads will cause internal stresses 
and strains to develop until an equilibrium point is reached. The result is an internal 
potential energy, or stored energy that is able to do work. In other words, potential strain 
energy is the potential energy due to internal stresses referenced from a zero stress state of 
a deformable body. For a system with no losses, potential energy equals the work put 
into the system. To develop the strain energy relationship, consider a general body with 
externally applied loads as shown in Figure 3.1 below. 



Given a traction force T acting on a portion of the body dS with a outward normal 
n and that displaces 8u, then the work increment is: 



To develop potential strain energy, the basic governing equations are shown 

below: 

Boundary Conditions: 




T 



n 



8u 



Figure 3.1. External Loads to a General Deformable Body. 




(3.1 ) 



5 



T x =o xx n x +a xy n y +o xz n z 
T y =o xy n x +o yy n y + o yz n z 



(3.2) 
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Equations of Equilibrium: 



I , 



dx 


+ 


dy 


d^xy 


1 


do yy 


dx 


r 


dy 


dG xz 


+ 


dc yz 



dx dy 



dz 

dz 

dz 



= 0 
= 0 
= 0 



Strain-Displacement Relationships (Engineering Strains): 



e 



xx 



£ 



•TV 



du 

dx 

dv du 
dx dy 



yy 



dv 

dy 

dw du 
dx dz 



yz 



dw 

dz 

dw dv 
dy dz 



(3.3) 



(3.4) 



B. DERIVATION OF STRAIN ENERGY RELATIONSHIP 



From this point on, tensor notation will be used as a shorthand notation to develop 
the strain energy terms. The governing equations are rewritten below. The only change 
is in the shear strain terms, in which the tensor shear strains are 1/2 of the value of the 
engineering shear strains. 

Boundary Conditions: T { = G ij n j ( 3.5 ) 



Equilibrium Equations: 



3a, , 

aT -0 '" = ° 



(3.6) 



Strain Displacement Relationships: 



e v = 2^ U ‘j +Uj ^ 



(3.7) 



Returning to equation 3.1, we substitute in the boundary conditions and expand: 

dW = JJ T, 8m, . dS 

S 

dW = jj(a t n J )Su l dS 

s 

S 
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At this point, the Gauss-Divergence theorem can be applied, and the internal 
strain energy can be obtained from the externally applied loads: 

= dv 

V 

V 

When the equilibrium equations are substituted into the above equation, the first 
term goes to zero. The result is expanded to form a symmetric and a skew-symmetric 
matrix: 

V 

SW = Jjj CT.J+fSu,, + &<,.,)+ ,(&<,., - Su (J )]dV 

5W = J|jo,(5e, J+ 5o), J )</r 



where to is defined as: 



to 



= \{ U U 



However, since to is a skew-symmetric matrix, and since a skew-symmetric 
matrix multiplied by a symmetric matrix is zero, the last terms disappear in the above 
equations. The final result is: 

SW = JJJ afadV (3.8) 

v 



C. STRAIN ENERGY DENSITY 

Equation 3.8 is the differential strain energy. To find the total strain energy, the 
integration along each strain can be accomplished inside the volume integral. Strain 
energy density is the strain energy per unit volume and can be represented as the 
integrand of equation 3.8. The 5 operator will be changed to a differential, and when 
considered from a zero stress state, this becomes: 

W.=j o ‘'<vfe, (3,9) 

where e,y is the strain value at the final stress state. 
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1. 



Strain Energy Density for the Elastic Case 



To calculate the strain energy density for the elastic case, the elastic stress-strain 
relationship, given below, is substituted into equation 3.9: 



a,=2Ge,+te iy 5, 



with X given as: X = ■ 



vE 



'• G= 2(T^) ; and 5 " = 



(l + v)(l-2v) 

The evaluation of the integral is straight forward: 

w:=i:(2cc M )^ 



[>:■' = j 

lO; i*j 






j, e • 

W e = Ge e + — — 

o V V 2 



5, 



W'/=y(2Ge, + Xe,,8„) 



(3.10) 



(3.11 ) 

(3.12) 

(3.13) 

(3.14) 



Replacing the term in brackets again by the elastic stress-strain relationship, the strain 
energy density can be written as: 

(3.15) 

Geometrically, this is the area under the linear stress-strain curve for the case of 
uniaxial tension. 



2. Strain Energy Density for the Plastic Case 

To evaluate the strain energy density for the plastic case, the strain energy density 
relationship, equation 3.9, is separated into an elastic term and plastic term: 

K = w; +wf=j o,*; +J < 3. 16 ) 

The first term is equation 3.15, and the second term can be manipulated in the same 
manner as Equation 2. 14 to give: 

W =—ae e +ae p -\e p dG r (3.17) 

o 2 tJ lJ */ v J v i j v ' 
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The stress-strain relationship for elasto-plastic behavior for deformation theory based on 
the Prandtl-Reuss equations is given as Equation 3.18 [Ref. 8]. 



£ .; = 



■ 1 


V 


3 1 




-Ha-5 ' 


— a 


a,, 5. 


_j_ 


G 


t_2G " 


g kk ij 


2 E, 


u 


3 V J 



Where E s is the secant modulus of the effective stress versus plastic strain curve: 



E =■% 



(3.18) 



(3.19) 



where the effective stress, a e , and effective plastic strain strain, e/, are defined as: 
Effective Stress: a e = ( 3.20 ) 



Effective Plastic Strain 



- P — ll c P e t 
=> ” V* e v e (i 



(3.21 ) 



where the stress deviator tensor is: s i: = a, ; - — —5:, 

l J ^ U 

Substituting the plastic strain component of the relationship given in Equation 3.18 into 
Equation 3.17, the strain energy density becomes: 



W o=\aijel+G u 



< 3 - 22 > 



If the Ramberg-Osgood uniaxial stress-strain curve is rewritten in terms of effective stress 
and effective strain, it becomes: 



£ = £ f +£ p =— + 
c f TC « £ 



K 



( 3.23 ) 



Substituting the plastic portion of the effective strain of the Ramberg-Osgood relationship 
in Equation 3.23 into the definition of the secant plastic modulus, Equation 3.19, it can be 
rewritten in terms of the effective stress and the Ramberg-Osgood material constants. 



E. 



V 



K) 



1 N \ n lltl 
a. n 



K 



(3.24) 



Substituting this back into equation 3.22 results in the strain energy density relationship 
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in terms of stresses only (based on the Ramsberg-Osgood stress-strain curve): 










2K" 





-—5 

3 * 



d<5„( 3.25) 



Before calculating the strain energy density from the finite element results, the 
finite element program as applied to notched geomentries will be verified in Chapter IV. 
Following verification of the FEM results, the strain energy density Equations 3.9, 3.17 
and 3.25 will be numerically integrated to calculate the actual strain energy density. 
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IV. VERIFICATION OF FINITE ELEMENT METHOD 



To evaluate the strain energy density concept for calculating stress concentrations 
in the plastic zone, a finite element model was used to obtain the stress and strain data as 
an input to the calculations. However, verification of the finite element data was 
paramount to ensure that a legitimate analysis is performed. The finite element program 
used was I-DEAS Master Series™ Release 1.3 and 2.1. Quadrilateral plane stress and 
plain strain elements were employed, and a non-linear stress-strain analysis was 
performed. 

The I-DEAS Master Series™ offers several approaches to createing a finite 
element model. The general approach was to create a three dimensional model to 
represent the physical specimen. Next, an element mesh was created by specifying the 
type of element to be used and the number of elements on each side of the surface being 
meshed. All plane stress and plane strain analyses throughout this thesis used an eight 
node quadrilateral element. These elements are two-dimensional, with nodal degrees of 
freedom consisting of translation in x and y directions, and rotation about the z axis. 

Only one face of the model was required to be meshed. A mesh refinement routine can 
be used to refine the element shapes to reduce distortion (i.e., skewing and stretching), 
and this routine was used to refine the finite element meshes for all analyses. After 
generation of the mesh, boundary conditions were applied to represent external loads and 
fixed displacements. Note that boundary conditions can be applied at lines of symmetry; 
thus reducing a block by one half or one quarter of the original geometry, resulting in 
significantly reduced computation time. 

Effective mesh generation is paramount in obtaining a correct solution for the 
finite element method. A coarse mesh (few elements, large distance between nodes) will 
not produce the correct solution while an overly fine mesh (many elements, small 
distance between nodes) may result in excessive computational time and storage. 
Additionally, a dense mesh near large stress gradients is not only necessary for accurate 
solutions, but also to provide the data for effective post processing analysis. Thus, for 
each problem, a mesh study was performed to determine the optimal mesh size and 
distribution to use. This involved running several finite element models with varying 
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mesh sizes and distributions, and determining at what point further mesh refinements did 
not effect the solution. 

Two additional factors effect the non-linear portion of the analysis: the non-linear 
stress-strain relationship and the iteration procedure. For a non-linear stress-strain 
analysis, the I-DEAS Master Series™ program uses input data points to model the stress- 
strain curve. The program assumes a constant stiffness between each point, breaking the 
curve up into individual linear sections. A maximum of twenty points can be entered. 

To ensure a smooth curve, all twenty points were used to define the stress-strain curve, 
with denser groupings at the highest curvature. Two basic iteration procedures are 
available: either the Newton-Raphson method or the Modified Newton-Raphson method. 
The Newton-Raphson method was used, and although this may mean more computational 
time per iteration, it generally iterates at a faster rate than the Modified Newton-Raphson 
method, requiring fewer iterations. For either method, a convergence requirement can be 
set, specifying the minimum value a function can change from one iteration to the next 
before the solution is considered to have converged. 

Two basic types of comparisons can be made: the first to an analytic solution, and 
the second to experimental data. The first type provides a comparison to an exact 
solution, and differences between the finite element method and an analytic solution 
should be minimal. The advantage of this type of comparison is that it will show the 
error of the finite element method compared to the governing equations of solid 
mechanics. The comparison of the finite element method to experimental data will 
indicate how well the finite element method models actual mechanical behavior. The 
critical area in this case is the elasto-plastic relationship, or how well the Prandtl-Reuss 
equations represent this relationship. This includes modeling the effective stress and 
effective strain with a uniaxial stress-strain curve. 

Full stress field measurements in solid mechanics must rely on surface techniques 
and assume plane stress or plane strain conditions. There are also approximations 
necessary in obtaining the elasto-plastic measurements, such as reducing the data via the 
Prandtl-Reuss equations. Another disadvantage of using experimental data is the percent 
error introduced when making measurements. Despite these difficulties and experimental 
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error, qualitative experimental results are indispensable, and when executed properly, can 
show were the mathematical model differs from the actual physics. Numerical 
computational techniques will have no validity unless they can be shown to truly 
represent the actual physical phenomenon. Therefore, despite the difficulty in obtaining 
experimental data, an attempt has been made to compare the finite element method to 
previously published experimental results. 

A. VERIFICATION OF ELASTIC FINITE ELEMENT MODEL 

1. Elliptical Hole in an Infinite Plate 

The first step in the finite element modeling was to ensure the program could 
obtain accurate results for a notched type of geometry under linear elastic loading. An 
elliptical hole in an infinite thin plate was used for this comparison. Durelli [Ref. 9] 
gives the stress distribution around the elliptical boundary and compares this with 
experimental data. Brown [Ref. 10] gives the complete stress field for any uniaxial load. 
The stress distributions for a infinite plate with an elliptical hole are given in elliptical 
coordinates, which are shown in Figure 4. 1 and Figure 4.2, for an ellipse of a = 1 and b = 
0.5. 




Figure 4.1. Elliptical Coordinate System. 
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£, = constant 




r) = constant 

Figure 4.2. Normal Stress Components on an Element Referenced to Elliptical 

Coordinates. 



The conversion from elliptical to Cartesian coordinates is given by: 

x = c cosh i; cos r| 
y = c sinh % sin r| 

where c is the focal distance of the ellipse given by: 

c 2 = a 1 —b 2 

The stress around the circumference of the elliptical boundary is: 

_ 21 sinh 2E„ + e 1(, ° cos 2ri - 1 
a = 5 e 2 ^ 2*— - '■ — 

n cosh 2l So - cos 2r\ 



(4.1) 



(4.2) 



were C,o is 
The stress 



the value of C, at the boundary of the ellipse. 

components, as shown by Brown [Ref. 10], throughout the body are given as: 



5 

0,1 ” 2 






cos 2r\ - cosh + cosh(2^ - 2^ 0 )cos 2r] 
cosh 2^ - cos 2x\ 



+ 



sinh 2^ 0 sinh 2% - sinh 2% cos 2r| 
(cosh 2^ - cos 2ri) 2 



(4.3) 



sinh 2^ - e^° 2 ^ sinh 2^(cosh 2^ 0 + 1) — sinh 2^ cos 2ri 



cosh 2% - cos 2t| 



(cosh - cos 2r|) 



and 
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(4.4) 



o f = — 






cos 2r\ + cosh 2\ n - cosh(2^ - )cos 2r\ 
cosh 21; - cos 2r| 



sinh 2% 0 sinh 21; - sinh 21; cos 2rj 
(cosh 21; - cos 2r|)~ 



sinh 21; - e^° 2 ^ sinh 2l;(cosh 21; 0 + 1) - sinh 2% cos 2r| 

cosh 21; - cos 2r) (cosh 2\ - cos 2r\f 



Along the x axis, = 0, r| = 0, therefore g x = and a y = On- On the y axis, a^n = 0, 
t) = n/2, therefore a y = and o x = 

To approximate an infinite plate, a plate 40 inches by 40 inches with a central 
elliptical hole with a major axis of 1 inch and a minor axis of 0.5 inches was used. A far- 
field load of S = 1000 psi was applied. Due to symmetry, only one quarter of the block 
was modeled (20 inches by 20 inches). The boundary conditions required to accomplish 
this symmetry were zero displacement in the y direction along the x axis, and zero 
displacement in the x direction along the y axis. This is shown below in Figure 4.3. A 



S 




Figure 4.3 Layout of model representing an elliptical hole in an infinite plate. 
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mesh of 30 by 30 quadrilateral plane stress elements was used. Figures 4.4 and 4.5 below 
show the layout of the finite element mesh used. 




Figure 4.4. Finite element mesh for elliptical hole in an infinite plate. 




Figure 4.5. Detail of finite element mesh around elliptical hole. 



Comparison of the finite element results with that of the analytical solution shows 
very good agreement between the two. Figure 4.6 shows the tangential stress along the 
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boundary of the ellipse plotted against the x-axis location. The finite element results 
coincide with the analytic solution across the complete boundary. 
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Figure 4.7 below shows the percent error from the analytic solution. The error is 
generally less than 1.0%, even at the point of maximum stress (x = 1). The greatest error 
occurs near the point of sign reversal (;t = 0.65) and is a maximum of 7.45%. However, 
the stresses at this point are on the order of 5 times less than the applied far-field stresses, 
and the actual solution crosses zero. Therefore, a better representation of the amount of 
error is to normalize the difference between the finite element data and the analytic 
solution by the far-field stress vice the analytic solution itself, as stated in Equation 4.5. 
This is shown in Figure 4.8. 

. . ,. , _ _ ^ Finite Element ^Amalytic Value 1A - „ . , ^ , 

Normalized Stress Error = x 100 % ( 4.5 ) 
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Figure 4.7 Percent Error of tangential stress <7 n along elliptical boundary. 




Figure 4.8 Percent Normalized Error of tangential stress o v along elliptical boundary. 

A comparison was also made of the stress distribution along r| = 0 and r\ = n/2. Figure 
4.9 below shows the Ct, stress distribution along r| = 0 for both the finite element 
calculations and the analytic solution. Very good correlation is shown, with the 
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maximum error at 1.23% (see Figure 4. 10). Figure 4. 1 1 below shows the Gr| stress 
distribution along T| = 7t/2 for both the finite element calculations and the analytic 
solution. Again, very good agreement of the finite element method with the analytic 
solution was achieved. The normalized error is shown in Figure 4. 12. 



i 




Figure 4.9. Axial stress ( O n ) along r] = 0. 
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Figure 4.10 Finite Element Error of Axial Stress ( (J^ ) Along rj = 0. 
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Figure 4.11 Axial Stress ( ) Along rj = n/2. 
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Figure 4.12. Finite Element Normalized Error of Axial Stress ( o n ) Along r\ = n '/2. 



2. Conclusion of Elastic Finite Element Model 

The comparison of the finite element results with the analytic solution of the 
stress distribution around an elliptical hole shows that the IDEAS finite element program 
can be used to accurately analyze notch stresses under linear elastic loading. Even though 
the maximum error along the hole edge was 7.45%, this occurred at a point were the 
stresses were only 3% of the far-field applied stress. The focus of these finite element 
analyses was to obtain accurate stresses and strains in the vicinity of the plastic zone. 
Therefore, this error at the low stress regions becomes insignificant when compared to the 
magnitude of the applied loads and of the peak stresses. 

B. VERIFICATION OF ELASTOPLASTIC FINITE ELEMENT MODEL 
1. Uniaxial Test of Stress-Strain Curve 

The first elastoplastic comparison of the I-DEAS™ Software was on a long narrow 
block to verify the input of the stress-strain curve. A finite element mesh of 3 elements 
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by 100 elements with the outer dimensions 1.0 m by 0.03 m, shown in Figure 4.13. This 
simple case illustrates that the elasto-plastic solution converges to the input stress-strain 
curve. Figure 4. 14 shows the stress-strain curve used, along with the FEM calculations. 



in mini 

Figure 4. 13. Finite Element Model of Uniaxial Tension. 




Figure 4.14. Results of Uniaxial Test Case. 



Infinite Plate with a Circular Hole 



a. Analytic Solution - Elasto-plastic Case 



The elastic solution of the stress distribution of a circular hole in an 
infinite plate with a far-field applied load of o r = o, as shown in Figure 4. 15, can be 
obtained in closed form and is shown in Equation 4.6. To obtain the elasto-plastic stesses 
and strains, both incremental and deformation theory have been applied to this problem 



0=0 



a 

r J 



o Q =o 






1 + - 



(4.6) 



26 



[Ref. 1 1, 12, 13, and 14]. Tuba [Ref. 12] and Chakrabarty [Ref. 14] show that both these 
solutions are identical. The solution below is taken from Chakrabarty [Ref. 14]. 




Figure 4.15 Infinite plate with circular hole 



Applying the governing equations in polar form for axisymmetric loading, 
one obtains the stress components in terms of constant c 2 and C 3 that can be solved for the 
imposed boundary conditions: 

a r =2 c 2 +y a e =2c 2 -^ (4.7) 



The elastic solution shows that at the hole edge, Gq is twice the far-field value a, with all 
other components being zero. From the von Mises yield criteria. Equation 4.8, yielding 
will occur once a is greater than one half of the yield stress Go , and a region of plastic 
deformation will extend from the edge of the hole out to a certain radius of r > a. If the 
radius of the elastic/plastic boundary is given as c then the boundary conditions for 
Equation 4.7 are G r = o rc at r = c and G r = G at r = °°, where o rc is the yet to be determined 
value of the radial stress at the elastic/plastic boundary. 

<V = <V-a 0 a r + <V (4.8) 
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The constant C 3 , in Equation 4.7 can now be solved for in terms of Go and c, and the 
stresses in the elastic region reduce to: 



g =o- 



Sr 1 



1- 



G a — G + 



GqC 



V3r 2 



1- 



(4.9) 



In the plastic region, the Von Mises stresses can be expressed in parametric form as: 



G r - ^ G e sin 4> 






G. COS 



K ^ 

■ T"* 

V6 



(4.10) 



The stress-strain constitutive relationships given in Equation 3.18 are restated in polar 
form: 



1 / \ 3 

i r =j{o r -vo 6 ) + - 



1 / \ 3 

ee=^K-vG r ) + - 



y E' s E 






(4.11) 



y E‘ s ~ E 






Se-^K+ae) 



The plastic secant modulus has been rewritten as follows: 



_L 

E. G„ 



e‘ - e 



El ~ E 



were Ej is the secant modulus of the total stress-strain curve. The solution of the stresses 
is independent of the elastic Poisson’s ratio, v, and therefore this can be set equal to 0.5. 
This transform simplifies the constitutive relations, while allowing solution for the 
correct stress values [Ref. 13]. With this simplification, and the substitution of the 
definition of the secant modulus as the ratio of the effective stress to the effective total 
strain, the stress-strain relations in the plastic region can be expressed as: 
k' 



£ r - £ < sin 



4>- 



£ 0 = £ COS(j) 



(4.12) 



where £ r and £ e are the strain values for an incompressible material (v = 0.5). Once the 



stresses have been found, then the actual strains can be found using Equation 4. 1 1 and the 
correct value for Poisson’s ratio v. 
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The strain £ e is related to the equivalent stress o e by Equation 3.23, or a similar 
stress-strain curve. Another such empirical relationship that can be used is a modified 
Ludwick curve [Ref. 14], which is separated into an elastic and a plastic range: 



o = 






a o 



E 

Or. 

£>— T 



(4.13) 



A plot of a family of curves of this form is shown below in Figure 4. 16, where Go = 



30,000, E = 30.0 x 10 6 , and n varies from 0 to 0.5. At the elastic/plastic boundary, the 




E 

Figure 4.16 Modified Ludwick Stress-Strain Curve. 



auxiliary angle 0 is found by setting the radial component of the stress in Equation 4.9 
equal with the radial stress component in Equation 4. 10. The substitution of r = c, and 
o e = o 0 results in: 



=sin~ 



So 1 



V 



2 a n 



--J1- 



(4.14) 



When stresses in Equation 4.10 are substituted into the equilibrium condition, the result 

is: 
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(4.15) 



1 do, 3d) 1 / <- x 

— — iukh-— =— (^- , an (i>) 

Likewise, when the strains in Equation 4.12, along with the stress-strain law of 
Equation 4.13, are substituted into the compatibility equation, the result is: 

1 do , 3d> n / r- \ 

— = n tan (b ~~~ + — (v3 tan 6 -3) 

G or dr 2r ' ' 



(4.16) 



The relationship between the auxiliary angle <|> , r, and n (material constant from 
Equation 4.13) can now be found by eliminating the equivalent stress from Equations 
4. 15 and 4. 16. The resulting differential equation is: 

3<|> (V3 - tan <)))(l + V3 n tan <f>) 

dr 2(l + n tan 2 <(>) ^ 



The solution to this can be solved by separation of variables, and can be reduced to r as 
function of 



—5- = -jj (V3 cos 0 - sin <}>)(cos <f> + V3rtCos(|)) 



-4 n 
l+3/i 2 



-V3- 



1 -n l 
' Mn 2 



(4.18) 



The equivalent stress must also be found, and this is obtained by eliminating the radius 
from Equations 4.15 and 4.16. The resulting differential equation leads to the solution of 



the equivalent stresses as a function of <J>, n, and a: 



1 do e ( V3 - tan <t> 
o e 9<(> ^1 + V3ntan<)> y 



(4.19) 



The relationship between 0 and o e is now found as: 



o e = K' 0 (cos 4> -t- -n/ 3 « sin <}>) 



-n(l+3n)/(t+3n 2 ) 



l+3n 2 



) 



(4.20) 



The stress concentration factor is obtained by setting <j) = <|> c and thus g((}>) = Go in 
Equation 4.20: 





V3(l-») 

2 



n(l+3n) 




V3 

e 






(4.21) 
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For a given applied load of o r = a, the plastic stresses in Equation 4. 10 can 
be obtained by first solving for (|) t . and K, as given in Equations 4.14 and 4.21. The 
auxiliary angle <|> in Equation 4.18 is solved for a given radius by an iterative method. 
Finally G e is solved from Equation 4.20, which is then substituted into Equation 4.10 to 
solve for the component stresses. For a work-hardening material ( n * 0 ), the strains are 
then found by Equation 4. 1 1. 

b. Analytic Calculations 

Two cases were solved analytically to be used with comparison to the 
finite element solution. The first involved an elastic perfectly plastic material that 
matched that used by Davis [Ref. 11] for his numerical results using incremental theory. 
The second case used the same elastic properties but changed the plastic property to a 
work hardening material. Table 4-1 below shows the material properties: 



Material Constant 


Value 




30,000 


E 


30 x 10 6 


V 


0.3 


G 


22,500 


n. Case 1 


0 


n, Case 2 


0.25 



Table 4-1 - Material Properties of Infinite Plate with hole 

For the two separate cases, the calculated values of <}) c , c, and K t are listed below. 
Figure 4.17 shows a plot of Equation 4.18 for the two cases. Note that there are two 
solutions shown, and care must be taken to ensure that the correct root is solved for. The 
mathematical software Maple V® was used to solve <J> for each input radius. This program 
allows the analyst to limit the range in which to search for a solution, and therefore only 
solutions for 0 < <|> < 60° were obtained without any difficulty. Once <]) was found for a 
specific radius, the calculations of G e , followed by G r and Gq » was straight forward. 
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Maple V 1 ® was programmed to make all these calculations at a radial distance 
corresponding to each node location along the x-axis. 



Variables 


Case 1 (n = 0) 


Case 2 ( n = 0.25 ) 


4>c 


0.32446 


0.32446 


C 


1.51550 


1.47886 


K t 


1.3333 


1.50195 



Table 4-2 Solution Variables for Case 1 and Case 2. 




Figure 4.17 Solution of Equation 4.18. 

c. Finite Element Solution 

The infinite plate with a circular hole was modeled using a 60 degree 
section on an annulus with a inner hole radius of a - 1.0 and an outer radius of 30. Due 
to the axisymmetric loading, a variety of finite elements are available that would 
efficiently model the axisymmetric problem; however, the above geometry with 
quadrilateral plane stress elements was used in order to verify these elements and the 
I-DEAS™ software for more general notched geometries. 

The first step was to determine beyond what mesh sizing the solution was 
mesh independent. Three separate meshes as listed in Table 4-3 were analyzed and are 
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shown in Figures 4. 18 and 4. 19. The coarse and medium element meshes each used 
equal numbers of elements on opposite sides of the geometry, while the fine mesh had 
one row of elements at the outermost radius that split into a finer meshing using 
triangular elements. This was possible since both the radial and cicumferential stresses 
rapidly approach the far-field value, and for the load used, reaches within 2% of the far- 
field value at a radial distance of 7.5a and within 0.5% of the far-field value at 15.13a. 
The finite element model extends radially to 30a, hence there is little gradient in the 
stresses and strains beyond a radial distance of 15a. 



Mesh 


Elements 


Nodes 


Coarse 


200 


661 


Normal 


800 


2,521 


Fine 


1378 


4,247 


Table 4-3. 


Mesh sizes for hole 


in infinite plate. 
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Figure 4.19. Exploded Views of Element Meshes Near Hole Boundary. 



The results of the mesh comparisons is shown in Figure 4.20. The mesh 
refinement results in only a slight change in the solution, particularly at the elastic/elasto- 
plastic boundary. The solution change from the medium to fine mesh is insignificant, and 
it can be concluded that the 20 by 40 element mesh provides a result that will not improve 
significantly with an increase in the mesh density. For the rest of the comparisons to the 
analytic solution, only the 20 by 40 element mesh will be used. 
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Figure 4.20 Circumferential stress for various mesh sizes. 



The results of the FEM were taken along the one line of elements at 9 = 0. 
As would be suspected, the angular change in the FEM solution was insignificant, as 
shown in Figure 4.21, which shows the effective stress distribution throughout the finite 
element model. The radius of the elastic/plastic boundary, c, corresponding to an 
effective stress of 30,000 psi, was interpolated between nodal points to be 1.5185 inches 
for case one, and 1.5 107 for case two, with the analytic values listed in Table 4-2. This is 
an error of only 0.2% for the former case, and 0.59% for the later. Note that the finite 
element value is linearly interpolated, and though the nodal values may be correct, this in 
itself will introduce a slight error. 

For both case one and case two, the finite element solution was found to 
match the analytic solution very well. Figure 4.22 shows the finite element stresses with 
the analytic curves superimposed for the perfectly plastic model (case 1). The normalized 
error was less than 1% for each data point (see Figure 4.23). The stress distributions for 
case 2 are shown in Figure 4.24. Again, very good agreement is made with the analytic 
solution, and the normalized error remains less than 1% for the work hardening material 
(see Figure 4.25). For the elastic/perfectly-plastic model, the strains in the plastic region 
cannot be obtained analytically. However, for case 2, one simply solves Equation 4. 1 1 . 
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The strain distribution for case 2 is shown in Figure 4.26. As in the case of the stresses, 
the finite element method produces very accurate results. The error was normalized with 
respect to ( S/E) as shown in Equation 4.22 for the same reason the stress error was 
normalized in the preceding section. 



Normalized Strain Error = 



C — p 

Finite Element Analytic Value 

(s/E) 



x 100% 



( 4.22 ) 



3. 38E+04 
3. 3 0E+04 
3. 20E + 04 
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2. 50E+04 
2. 40E+04 
2. 25E+04 




Figure 4.21. Effective Stress for Case 2. 
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Figure 4.22. Comparison of F 'EM stresses with analytic solution for case 1. 




Figure 4.23. Normalized, error of F EM stresses versus analytic solution for case 1. 
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Figure 4.24. Comparison of FEM stresses with analytic solution for case 2. 




Figure 4.25. Normalized error of FEM stresses versus analytic solution for case 2. 
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Figure 4.26. Comparison of FEM strains with analytic solution for case 2. 




Figure 4.27. Normalized error of FEM strains versus analytic solution for case 2. 



d. Conclusion on the Accuracy of the Finite Element Method 

Comparisons of the I-DEAS Master Series™ finite element program to the 
specific problems shown reveal the excellent results obtainable for both linear and non- 
linear analysis by the finite element method. The differences between the finite element 
calculations with exact analytic solutions for both stress and strain calculations were 
minimal, even in regions of high gradients. It was also shown that a mesh-independent 
solution was readily obtainable after one or two mesh refinements. From these results of 
the application of the I-DEAS Master Series™ program to model problems, follow-on 
finite element analyses for notch geometries were performed with a very high degree of 
confidence in the accuracy of the solutions. 

3. Comparison to Experimental Data 

a. Finite Element Modeling of Experimental Specimen 

Although the availability of recently published full-field stress data for 
basic notched and hole specimens under tensile loads is slim to non-existent, several 
works were completed more than thirty years ago using photoelastic techniques. Two 
notable works are from Durelli and Sciammarella [Ref. 15] and Theocaris and Marketos 
[Ref. 17]. A comparison to the Theocaris and Marketos experiment was chosen vice the 
Durelli and Sciammarella due to the fact that Theocaris and Marketos were able to show 
the maximum peak tensile stresses progressing away from the notch root. This peak o y 
stress progression occurs due to a multi-axial stress state and given stress distributions 
allowing for a higher o y stress before yielding occurs. This was also the result of all the 
finite element analyses completed. The test specimens used by Theocaris and Marketos 
were two sheets of aluminum alloy 57S, one with a hole diameter to width ratio of 1/2, 
the other 1/3. The later ratio of 1/3 was chosen, and the dimensions are shown below in 
Figure 4.28. 
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I 7 inches I 

f 177.8 mm 1 

Figure 4.28. Dimensions of Plate with a Central Hole, [After Ref. 17]. 



The stress-strain curve of the Aluminum alloy 57S as given by Theocaris 
and Marketos [Ref. 17] is shown in Figure 4.29, with a curve fit to the Ramberg-Osgood 
equation resulting in the coefficients listed in Table 4-4. 



Coefficient 


Value 


E (kg/mm 2 ) 


7,000 


K (kg/mm 2 ) 


33.1 


n 


0.048 



Table 4-4 Ramberg-Osgood Material Coefficients for Aluminum 57S. 
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£ ej f (mm/mm) 

Figure 4.29. Equivalent Stress-Strain Curve for Aluminum 57S, [After Ref. 17]. 

The first step, as in the case of the model of the hole in an infinite plate, 
was to determine at what mesh density mesh independence was achieved. Three separate 
meshes were generated and analyzed: a 20 by 20 mesh, a 30 by 30 mesh, and a 40 by 40 
mesh. These are shown in Figure 4.30. Figure 4.31 shows the y-component of stress for 
a load well into the plastic range. The change in the FEM stresses from the 30 by 30 
mesh to the 40 by 40 mesh are insignificant, hence the 30 by 30 mesh will be used for the 
comparison to the experimental data. 
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Figure 4.31. Distribution of G y for Various Mesh Sizes. 

Table 4-5 shows the applied loads used by Theocaris and Marketos. The 
first load resulted in the initial onset of plastic deformation, while the last load creates a 
plastic zone that extends a distance 3/4 of the hole radius from the edge of the hole and 
covering approximately 1/3 of the minimum cross section. 



Load Set Applied Force 
(N) 


S (MPa) 


Snet (MPa) 


I 


19162.97 


67.89 


101.84 


n 


21384.95 


75.76 


113.65 


ni 


25520.79 


90.42 


135.63 


IV 


30399.66 


107.70 


161.55 


V 


34286.34 


121.47 


182.21 


VI 


38926.48 


137.91 


206.87 



Table 4-5. Load history for Theocaris Photoelastic Experiment, [Ref. 17]. 
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b. Finite Element Results 

A comparison of the G y stress was made for all six load sets, and these are 
shown in Figures 4.32 and 4.33. Several differences are noticeable between the FEM 
results and the Theocaris and Marketos data. First, although Theocaris and Marketos 
show a peak G y stress that progresses inward as the plastic deformation increases, they 
also show that it initially decreases before reaching a maximum for the last three applied 
loads. The FEM shows that the maximum value of a y moves away from the hole edge, 
but it also shows that G y constantly increases until it peaks. The second difference 
between the FEM and the Theocaris and Marketos data is the magnitude of the decrease 
in c y near the edge opposite the hole. This results in significant disagreement between 
the two values, approaching an 80% difference for the first applied load. Figure 4.34 
shows the difference between the FEM calculations and the experimental data. One test 
of the accuracy of both results is to determine if equilibrium has been satisfied at the 
minimum cross section. The stress distribution curves were numerically integrated to 
determine the resulting force. These results are listed in Table 4-6. The values are only 
for half of the plate, hence the total applied force will be twice these values. Even though 
these calculations are only approximate, it is easily seen that the FEM has satisfied 
equilibrium, while the Theocaris and Marketos results have underestimated the stress 
distribution for the first three load sets. It should also be noted that the photoelastic 
analysis shows fringes in the regions of strain gradients, and since this is a region of 
relatively uniform stress and strain, the accuracy of the method degrades. 
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Figure 4.32. Distribution of c y for Theocaris model, load sets I to III. 




Figure 4.33. Distribution of <J y for Theocaris model, load sets IV to VI. 
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Figure 4.34. Difference of F EM and Theocaris data for a y . 



Load 

Set 


Applied 

Force 


FEM Calculations 


Theocaris Data 




F y at Far 
Edge (N) 


F y at minimum 
section (N) 


Error 


F y at minimum 
section (N) 


Error 


I 


19162.97 


19166.13 


0.02% 


17254.83 


-9.96% 


II 


21384.95 


21387.70 


0.01% 


20004.95 


-6.45% 


m 


25520.79 


25522.31 


0.01% 


24975.24 


-2.14% 


IV 


30399.66 


30401.50 


0.01% 


30473.38 


0.24% 


V 


34286.34 


34286.69 


0.00% 


35225.31 


2.74% 


VI 


38926.48 


38925.78 


0.00% 


39660.45 


1.89% 



Table 4-6. Equilibrium calculations at minimum cross-section. 



c. Conclusion on the Accuracy of the Finite Element Method 

The finite element analysis compared favorably with the experimental data 
of Theocaris and Marketos. Within 20 mm of the hole edge, the error was less than 10% 
for all but one load set. Additionally, the results at the edge itself were within 4%. In 
addition to the experimental errors referred to in the beginning of this Chapter, the data 
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used in this study was taken from published graphs, resulting in an additional error added 
to the comparison. Considering these shortcomings, the FEM analysis provided 
quantitative results as reasonably as could be expected, and matched the qualitative trends 
exceptionally well. 
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V. NUMERICAL CALCULATION OF STRAIN ENERGY DENSITY 



A. INTEGRATION ROUTINE 

In order to analyze the strain energy density from the finite element data, an 
integration routine must be used to calculate the plastic strain energy density as shown in 
equation 3.25. The accuracy of any integration routine will be dependent on the number 
of data points integrated and the scheme implemented. This is especially true at the 
‘knee’ of the stress-strain curve, or just past yielding where the curve bends. To increase 
accuracy, a third order accurate integration scheme was used by approximating the first 
and second derivative of each point by a central difference scheme, and using these values 
to perform the numerical integration. This scheme is similar to Simpson’s rule, but 
applied for non-equally spaced points. To develop the integration routine, a function I(x) 
is assumed in which: 




(5.1) 



Expanding I(x) about the / th point by Taylor series, we get: 




(5.2) 



were the subscripts refer to each data point as shown in Figure 5.1. 
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fix) 




Figure 5.1 Integration points for numerical scheme. 



However, from Equation 5.1, we know that: 

r = f(x) i"(x) = f\x) r\x) = r\x) 

Substituting these into Equation 5.2 and solving for the integration between the I th - 1 and 
the I th + 1 point and defining this value as I\ results in: 



f, = J f(x)dx 



x ,-[ 



I-, = I( X j + Ax,.)- I(x, - Ax,., )+H.O. T. 



/,= (Ax, + Ax, _,)/(*,) + 



A 2 * 2 
-^.-1 



2! 



/'(*,) + 



* 3 . A 3 \ 

Ax ( +Ax,._, 



/ 






3! 



f"(x t ) + H.O.T. (5.3) 



J 



The first and second derivatives of the function are approximated by a central difference 
scheme as follows: 



/'(*,) 



A x,-, 2 (/(*,-+, ) ~ f(x,)) + Ax ; 2 (/(x, ) - /(*,■_, )) 
Ax,-i AXj (Ax m + Ax, ) 



(5.4) 



and 



/"(*,) = 2 



A*i-i(/(*,+i ) ~ /(*,- )) + ) ~ /(*,- )) 

Ax^Ax^Ax,., + Ax,) 



(5.5) 
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Equations 5.4 and 5.5 are substituted into Equation 5.3, then summed over every other 
data point to obtain the final integration value: 



B. RESULTS OF NUMERICAL SCHEME ON UNIAXIAL CASE 

A comparisons of the above numerical integration routine to a standard 
trapezoidal integration routine were made, using the uniaxial strain energy density 
equation with the Ramberg-Osgood stress-strain relationship. The stress-strain curve 
integrated is shown below in Figure 5.2, with the limits for both stress and strain labeled. 
The material constants were n = 0.053, K = 95ksi, and E = 10 xlO 3 ksi. For the modified 
Simpson’s rule integration routine, two forms of the strain energy integral were used, 
Equation 2.6 and 2.14. The actual value is given by Equation 2.17. Note that for 
integration by the trapezoid rule, both of these forms of the strain energy density produce 
identical values; however, for the applied integration routine, an upper and lower bound 
is calculated. This is shown in Figure 5.3, which applied the routine to equally spaced 
points. Note that the actual implementation of this routine will produce improved results 
due to a more efficient distribution of the data points; i.e., few points along the linear 
portion and a closer distribution in the region of higher curvature. Since the modified 
Simpson’s integration routine produced better results than the trapezoidal rule for higher 
data points, this was chosen as the means of integration to calculate the strain energy 
density. Additionally, the basic form of the strain energy density integral as given in 
Equation 2.5 was able to be incorporated within the I-DEAS Master Series™ post- 
processing module; hence, this was the final form used to calculate the strain energy 
density. 




even i 



(5.6) 
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Calculated Stain Energy Density 
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Figure 5.2. Stress-Strain Curve for Integration Comparison. 
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Figure 5.3. Comparison of Integration Routines. 
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Error 



8 % 
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VI. CALCULATIONS AND RESULTS OF STRAIN ENERGY DENSITY 



In order to test the validity of the Glinka strain energy density proposal, the strain 
energy density for a given loading was calculated based on an elasto-plastic material ( W p ) 
and an elastic material (W e ), and comparisons between the two calculations were made. 
The strain energy density calculations were performed at each node throughout the FEM 
model, with the plastic strain energy density ( W p ) being numerically integrated at each 
load step as shown in Chapter V. In addition to the strain energy density calculations, 
comparisons were made between the notch stresses and strains based on the finite 
element method and values obtained via both the Glinka and the Neuber method. These 
comparisons will be made in Chapter VII. The finite element analysis was performed on 
a plate with symmetrical, semi-circular edge notches. 

A. NOTCH GEOMETRY AND MATERIAL SELECTION 

Two separate geometries were evaluated, one had a notch radius of 1.0 inch, and a 
plate width of 6.0 inches (r/D = 1/6), and the second had a notch radius of 1.0 inch and a 
plate width of 10.0 inches ( r/D = 1/10). These two layouts are shown below in Figure 

6.1, with the shaded regions representing the finite element geometry. For each 
geometry, plane stress (thin plate) and plane strain (thick plate) analyses were performed. 
The strain-hardening material used modeled 7075-T6 Aluminum. This was represented 
in the I-DEAS Master Series™ finite element program as 20 data points, shown in Figure 

6.2. This figure also shows the yield stress, determined to be 66 ksi based on a plastic 
strain offset of 0.002. However, it should be noted that the stress-strain curve departs 
from linearity at 40 ksi. The 20 data points were obtained by curve fitting the Ramberg- 
Osgood equation to actual stress-strain data from the Military Handbook V [Ref. 1 10]; 
then calculating these points from the resulting Ramberg-Osgood equation. Although this 
may have resulted in the stress-strain curve used by the FEM analysis differing from the 
actual 7075-T6 data, it ensured that the notch root strains and plastic strain energy density 
could be calculated from the finite element stresses using the Ramberg-Osgood equation. 
This enabled correlating predicted strains based on the Ramberg-Osgood equation and 
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either the Neuber or Glinka method to the finite element results. The material property 
coefficients, including values resulting from the data fit to the Ramberg-Osgood equation, 
are listed in Table 6-1. 
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Figure 6.1. Notch Geometry and Finite Element Boundary Conditions. 
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Coefficient 


Value 


E (ksi) 


10,000.0 


K (ksi) 


92.0 


Co (ksi) 


66.18 


n 


0.053 



Table 6-1. 7075-T6 Material Constants. 



B. FINITE ELEMENT MODELING 

Due to the two lines of symmetry for this geometry, the finite element model was 
reduced to one quarter of the physical model by applying appropriate constraints along 
each boundary, as shown in Figure 6.1. For each configuration, 21 increments were used 
to increase the loading from an initial nominal stress of 12 ksi to a final value of 49.5 ksi, 
or 75 percent of the yield stress. This load increment ensured convergence of the finite 
element solution and provided a small enough step to numerically integrate the strain 
energy density with reasonable accuracy. The first two load points were in the elastic 
range, and plastic deformation began at the third load step. The stress concentration factor 
K, was determined from the finite element analysis at the first load step, and found to be 
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2.042 and 2.421 for the narrow and the wide plate, respectively. Table 6-2 shows the 
load schedule for each of the geometries. 

As in the previous finite element analyses, a mesh comparison was performed for 
each of the geometries. The number of elements and nodes for each geometry used in the 
comparison are shown below in Table 6-3. For the wide plate ( r.D - 1:10), the medium 
mesh was obtained directly from the coarse mesh by an I-DEAS Master Series™ routine 
of splitting elements in the region of highest elastic strain energy density. The mesh 
layouts are shown in Figures 6.3 and 6.4. Comparisons of both stresses and strains were 
made for each of the mesh sizes at load step 21 . The comparison showed that the mesh 
refinement resulted in little change in the stresses and strains throughout the model. 
Figures 6.5 and 6.6 show the stress and strain distributions along the minimum cross 
section for the each plate. For the narrow plate at the notch root itself, the stress 
decreased by 0. 16% going from the coarse to medium mesh, and decreased by 0.03% 
going from the medium to fine mesh. For the strains, the medium mesh resulted in a 
0.24% decrease in strain compared to the coarse mesh, and there was no change 
compared to the fine mesh. Similar results are shown for the mesh refinement for the 
wide plate. The mesh adaptation used on the coarse mesh to obtain the medium mesh 
resulted in a decrease of the notch root stress by only 0.22%. The strain at the notch root 
decreased 0.58% from the medium mesh to the coarse mesh. These results showed that 
the medium mesh for both geometries was sufficient to provide mesh independent 
solutions. Hence, for all further finite element analysis, the medium meshes were used. 
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Load 


Notch r:D = 


1:6 


Notch r:D = 1 : 


: 10 


Step 


P (ksi) 


S (ksi) 


K? S (ksi) 


P (ksi) 


S (ksi) 


Kf (ksi) 


1 


8.00 


12.00 


24.50 


9.60 


12.00 


29.05 


2 


13.00 


19.50 


39.81 


15.60 


19.50 


47.21 


3 


15.00 


22.50 


45.94 


18.00 


22.50 


54.47 


4 


16.00 


24.00 


49.00 


19.20 


24.00 


58.10 


5 


17.00 


25.50 


52.06 


20.40 


25.50 


61.73 


6 


18.00 


27.00 


55.13 


21.60 


27.00 


65.36 


7 


19.00 


28.50 


58.19 


22.80 


28.50 


68.99 


8 


20.00 


30.00 


61.25 


24.00 


30.00 


72.63 


9 


21.00 


31.50 


64.31 


25.20 


31.50 


76.26 


10 


22.00 


33.00 


67.38 


26.40 


33.00 


79.89 


11 


23.00 


34.50 


70.44 


27.60 


34.50 


83.52 


12 


24.00 


36.00 


73.50 


28.80 


36.00 


87.15 


13 


25.00 


37.50 


76.56 


30.00 


37.50 


90.78 


14 


26.00 


39.00 


79.63 


31.20 


39.00 


94.41 


15 


27.00 


40.50 


82.69 


32.40 


40.50 


98.04 


16 


28.00 


42.00 


85.75 


33.60 


42.00 


101.68 


17 


29.00 


43.50 


88.81 


34.80 


43.50 


105.31 


18 


30.00 


45.00 


91.88 


36.00 


45.00 


108.94 


19 


31.00 


46.50 


94.94 


37.20 


46.50 


112.57 


20 


32.00 


48.00 


98.00 


38.40 


48.00 


116.20 


21 


33.00 


49.50 


101.06 


39.60 


49.50 


119.83 



Table 6-2. Load Increment Schedule. 



Geometry 


Mesh 


Element 


Nodes 


Narrow Plate 


Coarse Mesh 


298 


965 


(r:D= 1:6) 


Medium Mesh 


1,128 


3,507 




Fine Mesh 


2,291 


7,046 


Wide Plate 


Coarse Mesh 


456 


1,453 


(r.D = 1:10) 


Medium Mesh 


1,064 


3,233 



Table 6-3. Number of Elements and Nodes for Mesh Comparison. 
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Figure 6.3. Mesh Laly outs for Narrow Notched Plate. 
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igure 6.5. 



Distribution of (J y for Mesh Comparison of Narrow Plate FEM Model. 




Figure 6.6. Distribution of e y for Mesh Comparison of Narrow Plate FEM Model. 
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Figure 6.7. Distribution of o y for Mesh Comparison of Wide Plate FEM Model. 




Figure 6.8. Distribution of £ y for Mesh Comparison of Wide Plate FEM Model. 
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c. 



RESULTS OF FINITE ELEMENT ANALYSIS 



1. Computational Procedures 

The elastic strain energy density (VF e ) at each node was initially calculated at the 
first load step. To determine W e at later load steps, this value was simply multiplied by a 
factor of ( Pi/Po ) 2 , where P, is the load at the / th step, and P 0 is the initial load. Plastic 
strain energy density (W p ) was calculated incrementally with the integration procedure 
shown in Chapter V. To calculate a strain energy density increment A W p , two new data 
points were required, in addition to last point of the previous increment. Since the first 
load step was under elastic conditions, W p is equal to W e \ hence, the first load step 
resulted in the first plastic strain energy density calculation without requiring any 
integration. Subsequent calculations were performed at all odd load steps. 



2. Plane Stress Condition 



For both plate widths, the plastic strain energy density was found to be greater 
than the elastic strain energy density in the vicinity of the notch root. Figures 6.9 
and 6.10 are plots of W p at load step 21 for the narrow and wide plates, respectively, 
under plane stress conditions. These show that the plastic strain energy density 
throughout the model has its maximum value at the notch root, but rapidly approaches the 
far-field value away from the notch. Figures 6.11 and 6.12 show W p for several load steps 
along the minimum cross section of the plates, with W e included at the final load step. 

This shows that not only does W e give under estimated values at the notch root, but also 



W, 



ERROR 



w n -w 
w m 



( 6 . 1 ) 



follows a different distribution shape than W p . At the notch root itself for the final load 
step (S = 0.75Co), the difference between W p and W e reaches a value of 16.8% and 23.2% 
for the narrow plate and wide plate respectively. Figures 6. 13 and 6. 14 show contour 
plots of the difference in W p as compared to W e , as shown in Equation 6.1, so that positive 
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values indicate that the actual W p is higher than that estimated by the Glinka Method. 
These values are taken at the final load step (S = 0.75 Oo). In the plastic region, the 
greatest difference between the two calculations occurs not at the notch root itself, as may 
be expected, but slightly offset along the x-axis. This maximum error is also on the order 
of twice the value of that at the notch root. The location of the maximum difference 
between W p and W e corresponds fairly well with the location of the maximum o y value. It 
should be noted that although the contour plots show the maximum global error occurring 
slightly above the notch near the plate edge, this is also a region of very low to zero strain 
energy density, and therefore these errors are actually insignificant. Figures 6.15 and 6. 16 
show the error in strain energy density for each load step across the x-axis (at y = 0 ). 
These plots show that, as previously stated, the maximum error occurs at the notch root 
for small plastic yielding, then gradually progresses inward along the x-axis. Global 
errors for the previous load steps followed the save trends as that of the final load step. 
Starting with zero error at the initiation of plastic deformation, the regions of significant 
error (greater than 1%) start at the notch root, and move inward, while at the same time 
proceeding at an angle toward the vertical centerline of the plate. This trend also 
corresponds to the region of high C e and the growth of the plastic zone (this growth is 
well documented by Theocaris and Marketos [Ref. 18]). 
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Figure 6.9. Plastic Strain Energy Density W p for Narrow Plate in Plane Stress 

at Load Set 21. 



68 




Figure 6.10. Plastic Strain Energy Density Wpfor Wide Plate in Plane Stress at 

Load Set 21. 
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figure 6.11. Strain Energy Density along x-axis for Narrow Plate in Plane Stress. 
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Figure 6.12. Strain Energy Density along x-axis for Wide Plate in Plane Stress. 
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Figure 6. 13. Difference Between Wp and W e for Narrow Plate in Plane Stress at 

Load Step 21. 
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Figure 6.14. Difference Between W p and W e for Wide Plate in Plane Stress at 

Load Step 21. 
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FigUre 6 ‘ 15 ' Per ™™ Difference Between W p and W e for Narrow Plate in Plane 
Stress Along Minimum Cross Section . 
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Figure 6.16. Percent Difference Between W p and W e for Wide Plate in Plane 
Stress Along Minimum Cross Section. 
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3. 



Plane Strain Condition 



The results of the strain energy density calculations showed better agreement 
between W p and W e than for the plane stress condition. Figures 6.17 and 6. 18 are plots of 
W p at the final load step. As in the plane stress case, there appears a high peak value at 
the notch root, with a rapid drop-off to the far-field value. Figures 6. 19 and 6.20 show W p 
for several load steps along the minimum cross section of the plates, with W e included at 
the Final load step. Closer agreement between the two strain energy densities is shown 
than for the plane stress condition. The plane strain problem results in a configuration 
that is physically more constrained than that of the plane stress problem, and hence the 
amount of plastic growth at the notch under plane strain conditions will be less than that 
of the plane stress condition. Therefore, the basis that the strain energy density 
distribution in the plastic region remains relatively unchanged due to a high volume of 
elastic material surrounding the plastic region should be even more valid for the plane 
strain condition than that of the plane stress condition. This was in fact shown to be the 
case when comparing the finite element plastic strain energy density W p with that of W e , 
as can be seen in Figures 6.21 and 6.22, which are contours of the difference between the 
two energies. As in the case of the plane stress condition, W p for the plane strain 
condition is shown to be greater than that of W e . However, both the amount of error and 
region of significant error (greater than 1%), is much improved over the plane stress 
condition. At the notch root itself, the difference between the two calculations was only 
7.76% for the narrow plate and 10.6% for the wide plate. 
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Figure 6. 1 7. Plastic Strain Energy Density W p for Narrow Plate in Plane 

Strain at Load Set 21. 
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Figure 6.18. Plastic Strain Energy Density Wp for Wide Plate in Plane Strain 

at Load Set 21. 
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Figure 6.19. Strain Energy Density along x-axis for Narrow Plate in Plane Strain. 
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Figure 6.20. Strain Energy Density along x-axis for Wide Plate in Plane Stress. 
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Figure 6.21. Difference Between W p and W e for Narrow Plate in Plane Strain at 

Load Step 21. 
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Figure 6.22. Difference Between W p and W e for Wide Plate in Plane Strain at 

Load Step 21. 
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VII. STRESS AND STRAIN CALCULATIONS 



A. PLANE STRESS CONDITION 

1. Finite Element Method Results 

The stress distributions for plane stress for both geometries are shown in 
Figures 7. 1 through 7.4. These plots depict either G y or a x along the minimum cross 
section ( y = 0 ) as a function of the nominal loading. As stated in Chapter IV, the 
maximum a y value progresses inward from the notch root as yielding increases. This can 
be attributed to the O x component, which starts from zero at the notch root and rapidly 
approaches its maximum value inward from the notch root. This increase in the G x 
component results in a higher allowable a y than at the notch root before yielding occurs. 

If one considers the limiting case of a perfectly plastic material; then at the notch root, 
a y , would remain at the value of the yield stress after initial yielding. However, at some 
point inward from the notch root, a y would be higher than the yield stress due to the G x 
component. The other limiting case is an elastic only material, where a y would linearly 
increase with loading. As can be seen from the stress distributions, the results of the 
work-hardening material falls somewhere between these two extremes. 
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Stress, x-component (ksi) 




Figure 7.1. Distribution of C x Along Minimum Cross Section for Narrow Plate in 

Plane Stress. 
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Figure 7.2. Distribution of <3 y Along Minimum Cross Section for Narrow Plate in 

Plane Stress. 
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Figure 7.3. Distribution ofo x Along Minimum Cross Section for Wide Plate in 

Plane Stress. 
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Figure 7.4. Distribution of G y Along Minimum Cross Section for Wide Plate in 

Plane Stress. 
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2 . 



Notch Root Stress and Strain Calculations and Comparisons 



In addition to comparing the plastic strain energy density W p to the predicted 
strain energy density W e , comparisons were also made between the finite element stress 
and strain results with those predicted by the Glinka and Neuber methods. The form of 
the Glinka method used was Equation 2.18, while for the Neuber method, Equation 2.4 
was used. If the nominal stresses are high enough then S * Ee, and Equations 2. 18 and 
2.4 are not valid. They may be modified, however, by using the Ramberg-Osgood 
equation instead of Hooke’s law to determine the nominal strains [Ref. 1, page 138]. For 
example, the Neuber’s Method, from Equation 2.2 would become: 
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Likewise, the Glinka method, from Equation 2.12, would result in: 
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However, when these modified forms were applied to the loading levels analyzed in this 
thesis, the amount of improvement was minimal to none. For example, when this 
adjustment was calculated for the final load level in the case of the narrow plate under 
plane stress conditions, the error in notch root strains at the final load step decreased from 
- 1 1 .85% to - 1 1 .66% for the Glinka method, and for the Neuber method the error actually 
increased from 14.71% to 14.88%. 

In Chapter VI, the strain energy density calculations showed that the predicted 
strain energy density based on elastic material properties (W £ ) was less than the actual 
strain energy density ( W p ). From this comparison of strain energy densities, it was known 
that the Glinka method would under predict the stresses and strains at the notch root. The 
stress and strains were also predicted based on the Neuber method. As was stated in 
Chapter II, the Neuber method has been shown to overestimate the stresses and strains. 
This was also true for all of the configurations analyzed in this thesis. These results of 
the Glinka method under predicting the stresses and strains and the Neuber method over 
predicting the stresses and strains are shown in Figures 7.5 and 7.6. 
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7 igure 7.5. Stress/Strain Concentration Factors for Narrow Plate in Plane Stress. 
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Figure 7.6. Stress/Strain Concentration Factors for Narrow Plate in Plane Stress. 
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These plots show the stress and strain concentration factor as a function of the 
nominal load for the FEM analysis, the Glinka method and the Neuber method. 
Additionally, a stress based on the average of the Glinka and Neuber determined stresses, 
and a strain determined from the Ramberg-Osgood equation and this average stress was 
determined for all configurations. It should be noted that whenever strain values are 
based on an average of the Glinka and Neuber methods, this implies that it is based on the 
average of the calculated stresses, and not the average of the strains. The amount of error 
for the stress and strain predictions for both methods is shown in Figures 7.7 through 
7. 10. Note that the strain is more sensitive than the stress for both calculations. This is 
self-evident from the fact that for the uniaxial stress-strain curve beyond the yield point, 
strain is highly sensitive to changes in stress. This sensitivity is plotted in Figures 7.1 1 
and 7.12, which shows the percent error of the Glinka method predictions with respect to 
the percent difference in strain energy density. Note that this sensitivity is also dependent 
on loading condition as it relates to the current stress-strain relationship. For example, as 
the loading increases, the error in stress prediction appears to asymptote to a single value 
after an initial increase, while the error in strain is almost linear with respect to the error 
in strain energy density. As a material approaches perfectly plastic, then from Equation 
2.6, any change or error in W p will result in a linear change in £, since a will approach a 
constant value. 

For the notch geometry and plane stress condition, the Glinka and Neuber method 
give an upper and lower bound to both the stress and strain predictions. Results based on 
the average of the stresses of the Glinka and Neuber method are in good agreement to the 
FEM results. It should be noted after about a S/Go ratio of 0.5 to 0.6, the rate at which the 
error increases for the Glinka method appears to be constant, while even though the error 
for the Neuber method continually increases, that rate at which it increases diminishes. 
This results in the error for the average of these two methods to reach a maximum 
between a S/Go ratio of 0.54 to 0.64, then decrease as the loading increases. However, it 
does not appear to asymptote toward zero, but to merely change sign as the error from the 
Glinka method becomes greater than that from the Neuber method. 
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B. PLANE STRAIN CONDITION 



1. Finite Element Method Results 

The stress distribution for the plane strain case are shown in Figures 7.13 through 
7.18. These plots depict either a x , c y or a z along the minimum cross section ( y = 0 ) as a 
function of the nominal loading. When compared to the plane stress condition, the plane 
strain results in higher axial stresses. This is due to the o z stress, which results in a higher 
hydrostatic pressure for a given loading, thus reducing the amount of plastic growth. The 
reduced plastic growth likewise results in the stresses increasing at a higher rate than in 
the plane stress condition. Additionally, the plane strain condition results in higher stress 
gradients in the vicinity of the notch. As in the plane stress condition, the a y stress peaks 
at a point inward from the notch root as plastic growth occurs. For the plane strain 
condition, this also occurs for the ct z stress. 

2. Notch Root Stress and Strain Calculations 

For the plane strain analysis at the notch root, W p is a function of o y and e y . For 
the plane strain condition, however, £ y is a function of both a y and a z , and the strain 
energy density equation does not reduce to the simple uniaxial version. To solve the 
plane strain problem, Glinka [Ref. 4, 5, 6] uses the transformation as suggested by 
Dowling, et al. [Ref. 19] that relates the uniaxial stress-strain curve to a plane strain 
stress-strain curve. From this transformed stress-strain curve, £ y can be found directly 
from a y . To obtain £ y , £ y and a y are related to the uniaxial stress-strain curve as shown 
below: 
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The Ramgerd-Osgood relationship can then be fitted to these transformed stresses and 
strains, and new parameters are determined: 



£ >' = r + 



f <5 



K 



where E" = 






(7.4) 



For the 7075-T6 material parameters, the plane strain transformed parameters are 
shown below in Table 7-1. Once these values are found, the stresses at the notch root can 
then be found by substituting all the normal coefficients with the transformed coefficients 
in Equation 2.18 and Equation 2.4, then using Equation 7.4 to determine the strains. 



Coefficient 


Value 


E* (ksi) 


11,222 


K * (ksi) 


106.77 


n 


0.0541 



Table 7-1. Plane Strain Transformed Material Coefficients. 

Since the results of the strain energy density comparison was much improved for 
the plain strain condition, it was expected that the stress calculations using the Glinka 
method would also be improved. This was indeed the case, as is shown in Figures 7.19 
and 7.20„ which compares the stress concentration factor as a function of the nominal 
loading. Note that again an average value of the two methods was calculated, and plotted 
for comparison. Another plot of the strain values is shown in Figures 7.21 and 7.22. 

Here, KtS is plotted versus the calculated notch root strain £ y . Overlaid on this plot is the 
uniaxial stress-strain curve, the transformed plane strain stress-strain curve, and the notch 
stress versus notch strain results. Not only does this show that the Glinka method gives 
better results than the Neuber method, it also shows that the FEM notch root stresses and 
strains follow the transformed stress-strain curve. Therefore, one can conclude that the 
transformation used to obtain the plane strain results is valid at the notch root. The error 
in both the stress and strain calculations for both geometries is shown in Figures 7.23 
through 7.26. These show that the Glinka method provides results three times more 
accurate than the Neuber method. Additionally, strains calculated based on the average 
stresses of the two methods give slightly more accurate values than the Glinka method 
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itself. This process also results in strains that are slightly higher than those of the FEM 
analysis, vice those of the Glinka method, which are slightly lower than the FEM 
analysis. 
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Stress, x-component (ksi) 




Figure 7. 13. Distribution ofo x Along Minimum Cross Section for Narrow Plate 

in Plane Strain. 




Figure 7. 14. Distribution ofo y Along Minimum Cross Section for Narrow Plate 

in Plane Strain. 
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Stress, z-component (ksi) 




Figure 7.15. Distribution of <J , Along Minimum Cross Section for Narrow Plate in 

Plane Strain. 
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Figure 7.16. Distribution of O x Along Minimum Cross Section for Wide Plate in 

Plane Strain. 
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Figure 7.17. Distribution of a y Along Minimum Cross Section for Wide Plate 

Plane Strain. 
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7 igure 7.19. Stress/Strain Concentration Factors for Narrow Plate in Plane Strain. 




Figure 7.20. Stress/Strain Concentration Factors for Wide Plate in Plane Strain. 



96 



120 



100 



80 



go >> 

*£ *> 



60 - 



40 



20 - 



♦ o ▲ 
to ▲ 

♦ o ▲ 

>6 * 

♦o A 
♦O A 
^♦O A 
A 

€A 



,»^r 



•v 



ys 



S' 



oK 



O FEM 

♦ Glinka Method 

a Neuber Method 

7075-T6 Stress-Strain Curve 

— — Transformed Plane Strain Stress Strain Curve 

• Notch Root Stress-Strain 



-t — *- 



0.000 0.002 0.004 0.006 0.008 0.010 0.012 0.014 

8 (in/in) 



Figure 7.21. Stress and Stress Concentration versus Notch Root Strain E y for Narrow 

Plate in Plane Strain. 



C/5 
* 




£ (in/in) 



Figure 7.22. Stress and Stress Concentration versus Notch Root Strain e y for Wide Plate 

in Plane Strain. 
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gure 7.23. Error in Notch Stress (J y for Glinka and Neuber Method, Narrow Plate in 

Plane Strain. 
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Figure 7.24. Error in Notch Stress (Jyfor Glinka and Neuber Method, Wide Plate in 

Plane Strain. 



98 



15% n 



10 % - 



c 

<L> 

O 

Un 

l> 

CL, 



3 5% 



00 

o- 

o 

u- 

§ 

u 



0 % 



♦— Glinka method 
—* — Neuber method 

• — Average of Glinka and Neuber Stresses 



-5% - 
0.15 



• . • — »- 







0.25 



0.35 



0.45 

S/a 0 



0.55 



0.65 



♦ 

0.75 



Figure 7.25. Error in Notch Strain e y for Glinka and Neuber Method, Narrow Plate in 

Plane Strain. 




Figure 7.26. Error in Notch Strain £ y for Glinka and Neuber Method, Wide Plate in Plane 

Strain. 
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VIII. APPLICATION OF RESULTS TO FATIGUE CALCULATIONS 



A. STRAIN LIFE CALCULATIONS 

To determine the consequences of the Glinka method on fatigue life calculations, 
the results obtained in Chapter VII for the plane stress condition were applied to a strain 
life analysis as given in Equation 8.1 [Ref. 1]. From this relationship, the number of 
reversals to crack initiation ( Nf ) is based on the strain amplitude ( Ae ). The fatigue 
strength coefficient ( a/), fatigue ductility coefficient ( £/), fatigue strength exponent 
( b ), and fatigue ductility exponent ( c ) are material properties and are listed in Table 



9.1. For all fatigue calculations, it was assumed that the cyclic stress-strain curve 
remained constant (i.e., no hysteresis effects). Although this does not accurately depict 
the cyclic 7075-T6 aluminum properties, it does provide a means to compare the 
individual methods as they apply to fatigue calculations using the results shown in 
Chapter VII. Additionally, the fatigue life calculations are based on fully reversed 
loading, with the load levels as shown in Table 6-2 used as the amplitude of the 
alternating load. With these assumptions, the strain values obtained previously in 
Chapter VII were simply half of the strain amplitudes for the cyclic fatigue calculations. 




( 8 . 1 ) 



Coefficient 



Value 



c 



< 7 / (ksi) 

£f’ 

b 



191.0 

0.19 

-0.126 

-0.52 



Table 8-1. Fatigue Properties of 7075-T6 Aluminum. 
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B. RESULTS OF STRAIN LIFE FATIGUE PREDICTIONS 



Figures 8.1 and 8.2 show reversals to failure (crack initiation), Nf, as predicted by 
the FEM, Glinka method, and Neuber method. These figures show the predicted life as a 
function of the far-field loading. The fatigue life range covers from over 6,000,000 cycles 
at the initial load for the narrow plate, which can be considered an infinite life, to less 
than 1,000 cycles at the high loads. Figures 8.3 and 8.4 show the error in fatigue 
predictions based on the Glinka and Neuber results as compared to the FEM analysis. It 
is not surprising that the greatest error occurs at the higher loading, which corresponds to 
where the highest differences in strain calculations occur. Additionally, the Glinka 
method produces greater errors at the higher loads than the Neuber method. While the 
Neuber method appears to have a maximum error just under 40% between an applied 
load of ±30 ksi to±35 ksi, the Glinka method continually gets worse, obtaining a 91% 
error for the wide plate at fully reversed load of ±39.6 ksi. 
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Figure 8.1. Reversals to Failure for Narrow Plate in Plane Stress. 
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Figure 8.2. Reversals to Failure for Wide Plate in Plane Stress. 
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Figure 8.3. Error in Fatigue Life for Narrow Plate in Plane Stress. 
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Figure 8.4. Error in Fatigue Life for Wide Plate in Plane Stress. 

The comparisons of the Glinka method versus the Neuber method applied to 
strain life predictions show that although the error in strain calculations are of the same 
magnitude, underestimating the stresses and strains results in a greater error when 
calculating the fatigue life. Additionally, not only does this result in a greater error in 
fatigue calculations, but it errs on the high side of these life calculations. This could have 
alarming consequences when applied to safety critical parts. It should be noted that the 
strain-life calculations give only a model of when the actual crack initiation will occur. 
These comparisons do, however, provide an accurate picture of the trends in using either 
the Neuber or Glinka methods in a popular model that is used to make fatigue life 
predictions. 



104 



IX. CONCLUSIONS 



This study has examined the proposal by Glinka that the strain energy density at 
the notch root is the same regardless if the material is elastic or elasto-plastic. A detailed 
comparison of the two strain energy densities was performed not only at the notch root, 
but throughout the field of symmetrical, semi-circular double notched plates. These 
comparisons were made for both plane stress and plain strain conditions. Strain energy 
density was calculated based on finite element analyses that had been rigorously tested 
with analytic solutions and experimental data. The strain energy density was numerically 
integrated, applying 2 1 load steps to reach a nominal stress equal to three fourths of the 
yield stress for each configuration. 

In addition to strain energy density calculations, stress and strain calculations 
based on Glinka’s strain energy density proposal and the Neuber method were performed, 
and compared with the finite element method results. To calculate the notch root stresses 
and strains, a transformation to a plane strain stress-strain relationship was performed. 

For the plane stress condition, the strain results were applied to fatigue life predictions, 
using the relationship between number of reversals and the strain amplitude. 

A. FINITE ELEMENT CALCULATIONS 

Prior to calculating the strain energy density for the notched geometries, the finite 
element program I-DEAS Master Series™ was thoroughly tested by applying it to 
problems that offered an analytic solution and experimental results. The conclusion 
concerning the accuracy of the finite element program are listed below: 

• The finite element program provided very accurate results for the elastic 
analysis of an elliptical hole in an infinite plate under uniaxial tension. This 
included a comparison of all stress components along the axis of the plate and 
the edge of the ellipse. 

• The finite element analysis of a hole in an infinite plate under elasto-plastic 
loading again produced exceptional results when compared to the analytic 
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solution. This included the both stress and strain results in the plastic and 
elastic region. 

• The finite element analysis of the narrow plate with a central hole produced 
good results when compared to the experimental data of Theocaris and 
Marketos [Ref. 17], especially at the hole edge. When equilibrium at the 
minimum section was tested, the finite element analysis gave exceptional 
results, while the experimental data was in error by as much as 10%. 

The results show that the finite element analysis had been extensively verified 
with analytic solutions based on the governing equations of mechanics. The stresses and 
strains from the finite element program were input into an integration algorithm that 
calculated strain energy density. Since the stresses and strains have been shown to be 
highly accurate, it can be inferred that the strain energy density calculations are also 
accurate. 

B. STRAIN ENERGY DENSITY CALCULATIONS 

The strain energy density based on an elastic material and an elasto-plastic 
material was calculated at every second load step for a total of twenty one loads. A 
detailed comparisons of the two strain energy densities reveal that: 

• The strain energy density in the vicinity of the notch root based on elasto- 
plastic material properties is higher than the strain energy density assuming 
elastic only properties. 

• The plane strain condition results in better agreement between the two strain 
energy densities than the plane stress condition. 

• At the higher loads, the greatest deviation between the two energies occurs 
slightly inward from the notch root, vice at the notch root itself. 

• The distribution of the two strain energy densities along the minimum cross 
section not only differ in magnitude, but also in shape. 
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C. NOTCH ROOT STRESS AND STRAIN CALCULATIONS 

The notch root stress and strain values were calculated based on the Glinka 
proposal at the notch root, and compared to the finite element data and the Neuber 
method of estimating notch root stresses and strains. The comparison of these two 
methods revealed: 

• The Glinka strain energy density method under estimates the stresses and 
strains, while the Neuber method overestimates the stresses and strains. 

• For the plane stress condition, the two methods appeared to give an upper and 
lower bound. Taking the average of the stresses from each method and 
determining the strains from this value gave very good results. 

• The Glinka method, while under estimating the stress and strain values, 
produces results two to three times better than the Neuber method for plane 
strain conditions. 

D. IMPACT ON FATIGUE LIFE PREDICTIONS 

Fatigue life calculations were made based on the Glinka method, the Neuber 
method, and the average of the stress values of these two methods, and compared to 
fatigue life calculations based on the FEM results for the plane stress condition. The 
impact of the two methods is summarized below: 

• Since the Glinka method under estimated the stresses and strains, it 
overestimated the predicted life. However, the amount of error steadily 
increased with the applied load, growing up to 90% for the wide plate 
geometry at its final load level. For the Neuber method, the error appears to 
reach a maximum of 40% for a lower cyclic load and improves slightly as the 
load increases. 

• Estimating the fatigue life by using an average of the stresses of the two 
methods produce results that were ±10% from those based on the finite 
element results. 
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Even though the fatigue calculations were based on finite element data, with no 
comparisons made to experimental data, a valid comparison between the two methods 
was made. Since the Glinka method under estimates the stresses and strains, the fatigue 
life based on this method will be greater than the actual fatigue life. This would result in 
parts failing prior to their expected life cycle, and would be very detrimental unless safety 
factors were built into the structure and the calculations. 

E. RECOMMENDATIONS 

This study thoroughly analyzed the Glinka strain energy density proposal for 
stress concentration calculations at the notch root, and compared the findings with results 
from the Neuber method. For the geometries studied, it was shown that the Glinka 
method under-predicts the stresses and strains, while the Neuber method over-predicts the 
stresses and strains. Further comparisons should be made not only to different 
geometries, but also with different material properties. A study should also be performed 
by comparing actual fatigue data with predictions based on the two methods. Taking an 
average value of the two methods appears to give good results, and as in the case of either 
the Glinka or Neuber method, can be a means of quickly and efficiently computing notch 
root stresses and strains. This mean value of the two methods should also be compared 
with different geometries and material properties, along with comparisons to actual 
fatigue data. 
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